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1. INTRODUCTION 


With the use of composite (non-metallic) materials and 
microelectronics becoming more prevalent in the construction of both 
military arid commercial aircraft, the control systems have become more 
susceptible to damage or failure from electromagnetic transients. One 
source of such transients is the lightning discharge. In order to 
study the effects of the lightning discharge on the vital components 
of an aircraft, NASA Langley Research Center has undertaken a Storm 
Hazards Program in which a specially instrumented F106B jet aircraft 
is flown into active thunderstorms with the intention of being struck 
by lightning. The overall purpose of the program is to enhance the 
capability of detecting and avoiding the hazards associated with severe 
storms and improving design capabilities to protect aircraft systems 
from unavoidable hazards. 

One of the specific purposes of the program is to quantify the 
environmental conditions which are conducive to aircraft lightning 
strikes. To this end, the F106, through 1985, has made 1378 penetra- 
tions at altitudes ranging from 3.1 km up to 12.2 km (temperatures 
from +5°C to -55°C) and obtained at least 690 direct strikes to the 
airframe (Fisher, et al . , 1986). Most of the strikes have occurred at 
altitudes in excess of 6 km and temperatures colder than -20°C. Thus, 
a good data base has been established for the study of lightning 
interaction and environment at higher (colder) altitudes (tempera- 
tures). Analysis of these data has shown that the greatest risk of 
a lightning strike to an airplane occurs at altitudes between 11 and 
11.6 km (-40 to -45°C) where turbulence and precipitation intensities 
were classified as light (Fisher and Plumer, 1984). Analysis of UHF 
radar observations made during flight operations has indicated that 
these high altitude lightning strikes were triggered by the aircraft 
(Mazur et al . , 1984). 

The results from the Storm Hazards Program are in conflict with 
compilations of data concerning inadvertent lightning strikes to 
aircraft which show that most reported strikes occur in the vicinity 
of 4.5 km altitude, which has been very near the freezing level. 

Fisher and Plumer (1984) also state that the nature of the high alti- 
tude lightning strikes encountered (triggered) by the F106 have been 
generally of a low amplitude current nature, consistent with the known 
characteristics of intracloud lightning discharges (although some of 
the time-resolved characteristics appear to be different than antici- 
pated). They state the need for a larger data base for the analysis 
of lightning and strike characteristics at lower altitudes where these 
characteristics are anticipated to be somewhat different. The primary 
difficulty during the 6 field seasons of the F106 program has been the 
paucity of lightning strikes at altitudes below 6 km. To date, only 
75 direct strikes to the aircraft have been recorded below 6 km, 11% 
of the total, with 41 of these 75 occurring in the 1985 field season 
alone. Virtually all of these strikes occurred in the altitude range 



from 4.3 to 6 km (temperatures colder than 0°C) and were apparently 
triggered by the aircraft (a deduction based on the on-board sensors 
and UHF-band radar data). 

Given the need of obtaining a better data base on lower altitude 
lightning strikes, an effort was undertaken to determine possible 
methods of directing the F106 into lower altitude regions of thunder- 
storms where the probability of a lightning strike would be maximized. 
This effort involved the work of two organizations, the Institute of 
Atmospheric Sciences (IAS) at the South Dakota School of Mines and 
Technology and Electro Magnetic Applications, Inc. (EMA) of Denver, 
Colorado. This report summarizes the work accomplished by IAS 
personnel . 


2. OBJECTIVES 


The objectives of this research program involved the use of the 
IAS two-dimensional. Storm Electrification Model (SEM) to simulate the 
thunderstorm environment in which the F106 was operating. The results 
of the model simulations would be used in two ways. The first use of 
the model results was to provide EMA with data on the time evolution 
and spatial patterns of the electric fields, small ion concentrations, 
electric charges on cloud and precipitation particles, and the total 
charge patterns in the simulated cloud. EMA personnel would then use 
these data as initial and boundary values in their modeling studies of 
the aircraft lightning environment. 

The main use of the IAS model simulations was to analyze the 
results with an emphasis on looking for relationships between the 
electrical structure of the storm and the basic cloud structure. The 
intent of this analysis was to determine what observable character- 
istics of the cloud correlate with strong electric field regions at 
lower altitudes in the hope of devising a scheme for vectoring the 
F106 into such regions to increase the lightning strike probability. 

An additional objective of the research effort was to develop a 
lightning parameterization scheme for incorporation into the SEM. The 
SEM, in its original configuration, was only capable of simulating the 
development of the electrical aspects of a thunderstorm up to the time 
of first lightning. Without a scheme for simulating lightning, the 
buildup of the electric fields and charges would continue without 
limit, reaching unrealistic values and causing the simulation to 
terminate because of the necessity of reducing the model time step in 
order to maintain numerical stability of the electrical transport 
equations in these high electric fields. In nature, the electrical 
stresses that build up as the charge separation processes proceed and 
the electric field increases are relieved by the charge transfer that 
accompanies the lightning discharge. An analogue to the lightning 
discharge must be incorporated in the model in order to accomplish the 
charge transfer and relief of electric stress so that simulations may 
continue into the mature and decaying stages of thunderstorm evolution. 
Only in this fashion can the complete evolution of a thunderstorm be 
studied dynamically, microphysical ly, and electrically. 



3. MODEL DESCRIPTION 


3.1 The Base Model 

The theoretical framework is a deep-convection, slab-symmetric, 
two-dimensional, time-dependent (2DTD) cloud model which has been 
applied in the past to several atmospheric convective situations. 
Atmospheric wind, potential temperature, water vapor, cloud liquid and 
ice, rain, snow, and graupel/hail (in the form of ice pellets, frozen 
rain, graupel , and small hail) are the main dependent variables. The 
model has been developed from the works of Orville (1965), Liu and 
Orville (1969), Wisner et al . (1972) , Orville and Kopp (1977), and Lin 
et a! . (1983). A density-weighted stream function has been used to 
extend the model to deep convection. 

The nonlinear partial differential equations constituting the base 
model include the first and third equations of motion (from which a 
vorticity equation is derived), a thermodynamic equation, and water 
conservation equations (for the three phases). The model has been 
designed such that mesoscale convergence can be superimposed in the 
lower levels and divergence in the upper levels. The manner in which 
such convergence is applied to the model and further details of the 
hydrodynamic equations can be found in Chen and Orville (1980). 

The bulk-water parameterization used in this model divides water 
and ice hydrometeors into five classes: cloud water, cloud ice, rain, 
snow, and high density precipitating ice (graupel/hail) with exponen- 
tial size distributions hypothesized for the three precipitating 
classes. Figure 1 with an accompanying key in Table 1 shows the pri- 
mary cloud microphysical processes simulated in the model. Briefly, 
the production of rain from cloud water is simulated using equations 
based on the works of Kessler (1969) and Berry (1968). Graupel/hail 
Is generated by the aggregation of snow, by the capture of snow or 
cloud ice by raindrops (contact freezing), or by the probabilistic 
freezing of raindrops (Bigg, 1953) due to the inherent ice nuclei 
content of a specific volume of water at suitably cold temperatures. 

An approximation to the Bergeron-Findeisen process is used to trans- 
form some of the cloud water to snow. Growth of hail is governed by 
equations for wet and dry growth (Musil, 1970) and shedding of rain 
from hail is included. Cloud water may be transformed to cloud ice 
in the region between 0°C and -40°C using an equation developed by 
Saunders (1957). Natural cloud ice is normally initiated at tempera- 
tures of -20°C and colder, using an equation after Fletcher (1962) for 
the number of natural ice nuclei active. Homogeneous freezing occurs 
at -40°C. Accretional processes (including riming) involving the 
various forms of liquid and solid hydrometeors are simulated. 

Rain, snow, and graupel/hail have appreciable terminal fall 
velocities, while cloud water and cloud ice have zero terminal velocity 
and, hence, follow the airflow. Evaporation of all forms of hydrometeors 
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TABLE 1 


Synbol 

P IMLT 

P IDW 

P IHOM 

P IACR 

P RACI 

P RAUT 

P RACW 

P REVP 

P RACS 

P SACW 

P SACR 

P SACI 

P 8AUT 

P SPI» 

P SFI 

P SDEP 

P SSUB 

P SMLT 

P GAUT 

P GFR 

P GACW 

P GACI 

P GACR 

P GACS 

P GSUB 

P GMLT 

P GWET 


Key to Figure 1 


Meaning 

Melting of cloud ice to fora cloud water, T Tq. 

Depositlonsi growth of cloud ice at expense of cloud water. 

Homogeneous freezing of cloud water to fora cloud ice. 

Accretion of rain by cloud ice; produces snow or graupel depending on the amount 
of rain. 

Accretion of cloud ice by rqin; produces snow or graupel depending on the amount 
of rain. 

Autoconversion of cloud water to fora rain. * 

Accretion of cloud water by rain. 

Evaporation of rain. 

Accretion of snow by rain; produces graupel if rain or snow exceeds threshold 
and T < Tq. 

Accretion of cloud water by snow; produces snow if T < Tq or rain if T >_ Tq. 

Also enhances snow melting for T ^ Tq. 

Accretion of rain by snow. For T < Tq, produces graupel if rain or snow exceeds 
threshold; if not, produces snow. For T > Tq, the accreted water enhances snow 
melting. ~ 

Accretion of cloud ice by snow. 

Autoconversion (aggregation} of cloud ice to form snow. 

Bergeron process (deposition and riming) - transfer of cloud water to form snow. 
Transfer rate of cloud ice to snow through growth of Bergeron process embryos. 
Depositions! growth of snow. 

Sublimation of snow. 


Melting of snow to form rain, I ^ Tq, 

Autoconversion (aggregation) of snow to form graupel. 

Probabilistic freezing of rain to form graupel. 

Accretion of cloud water by gTaupel. 

Accretion of cloud ice by graupel. 

Accretion of rain by graupel. 

Accretion of snow by graupel. 

Sublimation of graiqpel. 

Melting of graupel to form rain, T > T n . (In this regime, P rArw is assumed to 
be shed as rain.) ~ 0 bACW 


Wet growth of graupel; may involve P^ and p gaci BUSt i»c ludc P GACW or 
P ^gp , or both. The amount of PQ^^WSich is notable to freeze is shed to 
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LlSLs. 1_: Cloud physics processes simulated in the model 
Table 1 for an explanation of the symbols. 
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and the melting of snow and graupel/hail are also simulated. For a 
more detailed discussion of this model, the reader is referred 
particularly to Orville and Kopp (1977) and Lin et al . (1983). 

3.2 The Electrical Model 

This base has been modified to include various processes involved 
in thundercloud electrification. These modifications involved the 
coupling of the electro-microphysical processes which account for the 
charge exchange between interacting particles within the cloud domain. 
The means of coupling the electro-microphysics to the dynamics and 
microphysics of the base model are based on ideas from Chiu (1978), 
Helsdon (1980), and Kuettner et al. (1981). Each of the five classes 
of hydrometeors is allowed to have a charge associated with it. This 
charge Is then transported with the species of hydrometeor according 
to the general equation 


^ *(§) inter • <» 


where Q represents the charge density (C nr 3 ) carried by the hydro- 
meteor class (positive or negative), v is the 20 velocity vector, p a 
is the air density, Ut is the mass-weighted mean terminal fall speed 
of the hydrometeor class, Kq, is the nonlinear eddy coefficient and 
(5Q/6t) inter represents the charge exchanged between classes of hydro- 
meteors due to interactions. The first term on the right is the 
advection term, the second is the fallout term (zero for cloud water 
and cloud ice charge), and the third is the eddy mixing term. The 
last term in (1) will be discussed below. 

The model framework accounts for the presence of small ions as 
well. The equation governing the number concentration of small ions 
is 


-|£i ,2 = -v.[n lj2 v ± ni s 2 Mi, 2 E - KmV n i,23 + 6 - an^ + Src - Sink. (2) 


Here n! 2 is the concentration of small ions (the subscript 1 repre- 
senting’positi ve ions and 2 the negative ionsi; y is the small ion 
mobility (m 2 V" 1 sec" 1 , pressure dependent); t is the 2D electric 
field vector; G is the ion generation rate due to cosmic ray ioniza- 
tion (ions m" 3 sec" 1 , dependent on height); a is the ion-ion recom- 
bination coefficient (1.6 x 10" 12 m 3 sec" 1 ); Src is the source term 
for small ions from processes other than cosmic ray ionization; and 
Sink is the sink term for small ions from processes other than recom- 
bination. An example of a source process for small iAns would be the 
evaporation of a charged hydrometeor, while a sink process would be 
Wilson capture (ion attachment to hydrometeors) . 



With the various types of charges accounted for in the model, a 
total charge density can be diagnosed at each grid point from 


p = e(nj-n 2 ) + EQ 


(3) 


where p is the total charge density (C m“ 3 ) and e is the electronic 
charge (1.6 x 10" 19 C). The first term represents the net charge due 
to the difference in small ion concentrations and the last term repre- 
sents the net charge on the five classes of hydrometeors. With the 
net charge determined, we then use Poisson's equation 


tf 2 <j> - p/ e 


(4) 


to determine the scalar electric potential, <f> (Volts). In (4), e is 
the electrical permittivity of air (8.86 x ID -12 kg -1 m~ 2 sec 2 C 2 ). 
Finally, the 2D electric field vector can be determined from the 
potential via Gauss's law 


E = 


-V<t> 


(5) 


The last term in (1) represents the charge exchanged between any 
two classes of hydrometeors during an interaction between those two 
classes. Such interactions may be accretional or non-accretional in 
nature. It should be noted that, in such interactions as P-jacr in 
Fig. 1 (Table 1), where cloud ice is nucleating rain to form snow, 
the charge on the cloud ice and rain are summed and the total is 
introduced into the snow charge resulting in a three bodied type of 
transfer. 

In the collisional process where a finite separation probability 
exists, two types of electrical interactions are possible: inductive 
processes, where the ambient electric field strength and its orien- 
tation with respect to the particles influence the magnitude and sign 
of the charge transfer; and non-inductive processes, where the ambient 
field exerts no influence, but such factors as temperature, liquid 
water content, and impact speed determine the charge transfer. The 
model has been configured to allow for the accounting of both induc- 
tive and non-inductive processes either individually or in concert. 
This is accomplished by starting with the following representation 
of the charge transferred between two interacting particles, 

+ 2 

Aq = AQ + 4ney t |E |cos(E,r|_s)r s + Aqi_ - Bq s , 


( 6 ) 



where AQ is the charge transferred due to any non-i nducti ve process, 

Yi is a dimensionless variable which depends on the ratio of the radii 
of the two interacting particles, T*ls is the vector along the line 
joining the centers of the two particles at the time of impact, r s is 
the radius of the smaller particle, qL (q s ) is the charge already 
existing on the large (small) particle and A (=1-B) is a dimensionless 
variable which depends on the radius ratio of the two hydrometeors. 

The first term, then, represents the non-inductive charge transfer, 
while the second term represents the inductive transfer. The last two 
terms account for the effect of charge already residing on the hydro- 
meteors. Switches are used to activate/deactivate either process. 

A single large hydrometeor may interact with many smaller 
particles in a unit time, so we Integrate over the volume swept out by 
the large particle 


ll l = jEf|V S L|N s AqS(a)dA (7) 


to obtain the charge acquired per unit time by the large particle. In 
(7), Ef is the collision efficiency between the two particles, ( Vjl | 
is the relative impact speed, N s is the number concentration of the 
smaller particles, Aq comes from (6), and $(a) is the angle-dependent 
separation probability. The integration results in three possible 
charging rates for the larger hydrometeor depending on which charge 
transfer mechanisms are active 


{ qi_ = 

At 


{c, + Bq s -AqJ| 


(8) 


where Ci = Ef |*^SL |N s llrL<S> 

^ 1 1 » A /C\ 4 


(ri_ being the radius of the large par- 


dcpa i at iuii 


k 4 1 \ 

pi uuci u i i i ujr ) 


c 2 = 


-AQ 


4lleYi |E |cos(E,VsL)rs<cosa> 


non-i nducti ve, 
inductive. 


1^411 e y i |E |cos(E,VsL) r s <COSa>_A Q combined. 


(9) 


Note that in (9) there are two important angles influencing the charge 
exchange: 1) the angle between the electric field vector and relative 
velocity vector of the two particles [cos(E,vsl)] and 2) the average 
coll i sion opnt act angle (<cosa>). For all practical purposes, the 
quantity |E| cos(E,v$l) is adequately approximated by -E z , the 
negative of the vertical electric field component. 


Since (8) implies that multiple charging interactions are possible 
for a precipitating hydrometeor per unit time of fall, and since the 
charge on the large particle appears on the right-hand side of the 
equation, it is not safe, a priori , to assume that 6qj_/6t is constant 
over a given time step. To account for this, we integrate (8) over 
the length of a time step. At, to obtain the charge transferred to the 
larger particle assuming that all quantities on the right-hand side 
are constant except q|_. This mav seem like a contradiction at first 
since both the electric field, | t | , and the charge on the small par- 
ticles, q s , appear on the right-hand side and would seem to be subject 
to time variation as well. It is safe to assume the constancy of q s , 
since it is highly unlikely that any one cloud droplet or ice crystal 
will undergo more than one non-accretional collision per time step. 

In the case of the electric field strength, it is a much more macro- 
scopic quantity, being determined by the entire charged volume of the 
cloud and not so significantly by the local charge changes exhibited 
by individual particles, although the effects of local volumes domi- 
nate. In spite of this, we do admit that the field can change rapidly 
on a local scale, but one other factor exists to help substantiate our 
locally steady assumption. This is the fact that the time step in the 
model is dependent on the electric field strength through the field's 
interaction with the small ion flux term. In order to maintain the 
stability of the ion transport equations, the time step must be 
reduced as the field strength increases. By the time the electric 
field has reached significant levels, the time step has been reduced 
to values of one second or less. Thus, we apply the steady field 
assumption with some confidence. 

Carrying out the integration of (8) results in the following 
expression for the charge accumulated by a precipitating hydrometeor 
during a time step At, 


- (Om-qLo) (l-e‘ At/Tl ) 


( 10 ) 


where q|_o is the initial charge on the large particle, and the 
parameters and tj are given by 


(-AQ + Bq s )/A 

-»• + 2 

Qm= < (4neyicos(E ,V|_s)r s <cosa> + Bq s )/A 


and 




(4Il£Y 1 cos(E,V$L)r s <cosa> - AQ + Bqs)// 1 


T i = (Ef |V S L |N s TT r L <S>A) 


-l 


non-inductive, 

inductive, 

combined , 


(ID 


( 12 ) 
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Finally, the charge accumulated by one precipitating hydrometeor 
is multiplied by Nl, the number concentration of the large particles, 
in order to arrive at the total charge density created by interac- 
tions, which appears as part of the last term in (1). The other 
components of the last term in (1) include the accretional terms, 
evaporational processes, and ion attachment processes. Examples of 
the formulation of the attachment and accretional processes are 
available in Chiu (1978) and Helsdon (1980). 

3.2.1 Charging processes 

A great deal of laboratory and theoretical work has been, and 
continues to be done to investigate possible mechanisms of charge 
separation. A non-exhaustive list of possible mechanisms includes: 

1) selective ion capture (Wilson, 1929); 2) convective charging 
(Vonnegut, 1955); 3) induction charging involving non-accretional 
collisions between various classes of particles in the electric field 
(Sartor, 1961, 1967, 1981; Muller-Hillebrand, 1954; Paluch and Sartor, 
1973); 4) effects of freezing potentials during wet growth of hail or 
splashing collisions between graupel and raindrops (Workman and 
Reynolds, 1948, 1950; Latham and Warwicker, 1980; Shewchuk and 
Iribarne, 1971); 5) non-inductive charging during the riming of 
graupel (Reynolds et a!, , 1957; Buser and Aufdermaur, 1977; Takahashi, 
1978; Gaskell and TTlingworth, 1980; Jayaratne et al . , 1983); 

6) riming of graupel followed by splintering (Latham and Mason, 1961; 
Hallett and Saunders, 1979); and 7) evaporative charging of ice 
crystals in penetrative downdrafts (Telford and Wagner, 1979). It is 
entirely possible that all of these mechanisms (or some as yet 
undiscovered mechanism) are active to one degree or another in various 

thunderstorm situations. However, in order to minimize the degree of 

complexity involved in the modeling work, it was necessary to limit 
the mechanisms included in the model to a workable set. 

The basic charge separation mechanisms included within the model 
can be broken down into four basic categories: 

1) Graupel interacting with particles in a dry growth mode. 

2) Graupel interacting with particles in a wet growth mode. 

3) Rain interacting with cloud droplets. 

4) Small ions interacting with all classes of hydrometeors. 

Within these various categories, the two types of electrical inter- 
actions discussed above may be possible, i.e., inductive and/or 
non-inductive transfers. 

The inductive water-water charge transfer interaction is generally 
believed to be ineffective in producing electrification typical of 
thunderstorm conditions on its own. Jennings (1975) has shown that 
such a process is capable of producing early electrification in clouds 
with fields up to 30 kV/m being possible. After this threshold is 



reached, the mechanism extinguishes itself due to the presence of such 
a field causing all collisions to result in permanent coalescence. 

This ability to generate early electrification well below thunderstorm 
values, but significantly above the background fair weather field, may 
affect the nature of other electrification mechanisms which must await 
the formation of ice particles for their initiation. 

The charge separation resulting from the interaction of liquid 
drops and ice particles may be of a non-inductive or inductive nature. 
Kuettner et al . (1981) have summarized the controversy regarding the 
action of a non-inductive, ice-water charging mechanism. They chose 
to include such a process in their model (using conservative values 
for the charge transferred per collision) based on the lack of conclu- 
sive evidence for eliminating the mechanism and the fact that the 
riming process is a well established microphysical phenomenon in 
thunderclouds. 

The non-inductive, ice-ice process is the subject of as much 
controversy as the water-ice mechanism; however, the controversy does 
not center around whether the process is effective, but rather exactly 
what physical mechanism is responsible for the charge exchange and, to 
a lesser degree, the magnitude of the charge transferred. Results 
have been summarized and assessed by Gross (1982). 

The inductive mechanisms involving water-ice and ice-ice 
interactions have been argued as being ineffective in nature because 
only grazing collisions near the electrical equator result in separa- 
tion of particles and charge in the ice-water case, and because of a 
relaxation time limitation on charge migration in the ice-ice case. 

In the ice-water case, it is true that grazing collisions will result 
in minimal charge separation; however. Sartor (1981) has shown that, 
due to the surface roughness of graupel particles, cloud droplets can 
bounce off of such rough particles from almost any impact angle 
negating the above restriction. In the ice-ice case, the relaxation 
factor may be included In the formulation, and the mechanism's impor- 
tance may be tested without great difficulty. In addition, the 
results of Tzur and Levin (1981) and Kuettner et al . (1981) have 

indicated that the combination of inductive and non-inductive 
mechanisms involving the ice phase yield the most realistic results. 

With this discussion in mind and referring to Fig. 1 and Table 1 
for the interactions listed below, the model in its present con- 
figuration can account for both inductive and non-inductive charge 
transfers when rain (Pgacr)» cloud ice (Pgaci)* and snow (^gacs) 
interact with graupel in category (1). In category (2), graupel can 
interact with rain and cloud water (P yacr and Pgacw through P ywe t) 
with the result that excess water is shed as rain (the wet growth 
process). In this category, only inductive charge transfers are 
accounted for due to uncertainty about the actual microphysical 
interactions taking place. This category is under further develop- 
ment at this time with a formulation for a non-inductive transfer 
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being devised. In category (3), a limited inductive mode transfer is 
allowed. In addition, both inductive and non-inductive charge trans- 
fers are allowed when cloud ice interacts with snow. The category (4) 
interactions involving small ions have been described above. 

3.3 Initialization 

3.3.1 Environmental initialization 

The base state of the atmosphere for this model is taken from a 
rawinsonde sounding typical of the type of day being studied. This 
sounding of temperature and moisture is, hopefully, representative of 
the atmosphere near the time of the formation of an observed cloud 
which may be compared with the model simulation, if such observations 
are available. The only modification made to the original sounding 
is to make the lowest layer adiabatic (if it is not already in such a 
state). The atmospheric winds are determined by projecting the winds 
from the sounding in the direction of the storm motion, since the model 
is two-dimensional. In addition, the winds determined by this pro- 
cedure are frequently reduced to some percentage of their original 
value because the use of full winds has a tendency to propagate the 
simulated cloud off the grid too rapidly. Dynamically, this approach 
is not on firm ground, and it would be better to subtract a mean wind 
from the actual winds and thus allow the model domain to translate. 
This approach will be incorporated in future work. 

In order to initiate convection, we use a combination of random 
perturbations in temperature and water vapor in the lower 3 km of the 
grid (with maximum amplitude ±0.5°C and ±7%) and a warm bubble in the 
vicinity of the domain center. The bubble has a maximum deviation of 
1.5°C, is 4.8 km wide, and is present between 400 m and 2 km height. 

3.3.2 Electrical initialization 

We assume that an initial steady state exists with ion production 
due to cosmic ray generation, ion loss due to recombination, and ion 
transport due to conduction in the ambient electric field all 
balancing. We ignore ionic diffusion and invoke horizontal homo- 
geneity so that the initial values of the variables are functions of 
height only. Then the aforementioned steady state balance can be 
described by 


•gy (u^^z) = G(z) - an t n 2 (13a) 


and 
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(13b) 


-d 

<JI 


(y 2 n 2 Ez) _ G(z) ~ an^ 


The terms in these equations have been previously defined. E z is 
the assumed fair weather electric field profile patterned after Gish 
(1944) which exhibits an exponential decrease with height such that 


E z = E 0 [b iexp(-a iz) + b 2 exp(-a 2 z) + b 3 exp(-a 3 z)] . (14) 


In (14), the parameters E 0 , b-j , and ai are varied to adjust the ver- 
tical profile. This vertical electric field profile can be related to 
the small ion densities by 


dE z 


dz 


e(n 1 -n 2 ) 

e 


(15) 


We subtract Eq. (13b) from (13a) and integrate in the vertical to 
obtain 


e(uinj + u z n 2 ) “ *T^z " ^c 


(16) 


where Xj is the total conductivity of the air and J c (A m -2 ) is the 
fair weather air-earth conduction current, which is assumed constant 
with height. Equations (15) and (16) can be solved simultaneously for 
ni and n z yielding 


n l 



1 

e(wi + t'2) 


(17a) 


and 



(17b) 


If we assume, following Shreve (1970), that the polar ionic mobilities 
vary exponentially with height such that 
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Ml, 2 = c 1>2 exp(f?z) 


(18) 


* 


where Cj = 1.4 x 10“ 4 and c 2 = 1.9 x 111 -4 m 2 V -1 sec -1 , and 
M 1.4 x 10 -4 m -1 , then (17a) and (17b) define the initial profiles 
for positive and negative small ions. Finally, either (13a) or (13b) 
may be solved for the galactic cosmic ray generation function G(z) 
which, after some algebraic manipulation, yields 


G(z) = eT^7Tn) t - 0EzEz + ^ Ez ^ 2 + EzEz] + anin2 


(19) 


The initial fair weather electrical state of the atmosphere is 
determined from Eqs. (14) and (17) - (19) by choosing appropriate 
values of E 0 , ai , b^ , and J c . 


3.4 Boundary Conditions 


The top boundary is assumed to be rigid with all variables held 
constant. The vorticity, vertical velocity, rain, snow, graupel , cloud 
water, and cloud ice as well as their associated charges are all set to 
zero. The stream function, entropy, water vapor mixing ratio, small 
ion concentration, and electric potential are maintained undisturbed 
at their initial values. The potential at the top of the model is 
obtained by integrating (5) over the depth of the domain and employing 
(14) for the electric field profile. This constraint on the electric 
potential is somewhat unrealistic and a boundary condition involving 
the derivative of the potential at the top is being devised. 


At the lower boundary, the vertical velocity, vorticity, and 
stream function are set to zero. The electric potential is also set 
to zero, consistent with the assumption that the earth's surface is a 
conducting plane. Evaporation and heating rates at the surface are 
prescribed as functions of time. Heat and water vapor are allowed to 
diffuse into the lower boundary. Clouds are not permitted to form at 
the surface, but precipitation is allowed to fall through the surface 
level. Cloud shadow effects (on heating) are also simulated at the 
lower boundary. 


At the lateral boundaries, the horizontal gradients of the stream 
function and the potential are assumed to be zero. Diffusional 
transport also assumes horizontal gradients of zero for all variables 
at the lateral boundaries. Both inflow and outflow from the domain are 
allowed. The constraints on the potential at the various boundaries 
result in the electric field components E x = 0 at all boundaries and 
E z being calculated at all boundaries. 


3.5 Numerical Techniques 


The equations are solved over a 19.2 km x 19.2 km domain with a 
200 m grid interval in both the X and Z directions. The advection 
technique used is that of Crowley (1968), which is first-order accurate 
in time and second-order in space. A two-step advection scheme is used 
(Leith, 1965) with vertical advection calculated first, followed by 
horizontal advection. Model variables are held constant at lateral 
inflow boundaries, whereas extrapolation via upstream differencing is 
employed at outflow boundaries. The general purpose Helmholtz solver 
in Cartesian coordinates (from the NCAR program library) is used to 
solve the Poisson-type equations for stream function and electrical 
potential. Centered differences are used throughout except at the 
upper and lower boundaries where second order, one-sided differences 
are used. 


4. RESEARCH RESULTS 


4.1 Wallops Island Simulation 

The first use of the SEM in connection with the Storm Hazards 
project was to simulate the 12 September 1983 operational flight of the 
F106 (flight 83-053) which experienced two lightning events near 28 kft 
during separate cloud penetrations. This simulation is termed a blind 
case because there were no specific observations with which the model 
results could be compared for verification. 

The SEM was initiated using the Wallops Island, Virginia, 00Z 
sounding from 13 September 1983, shown in Fig. 2, which was near the 
time of the flight operations (2215-2258 GMT) and was considered 
characteristic of the atmospheric state in which the storms were 
developing. The temperature has been modified in the lowest hundred 
millibars to make the profile adiabatic, as described in Section 3.3.1. 
Examination of the sounding shows that the Lifted Index is -6.7 and 
the K Index is 42.1, indicating a very high potential for strong 
thunderstorms on that day. In fact, severe thunderstorms might have 
been expected except for the fact that the upper level winds were 
quite light and exhibited little shear. Much convective activity was 
actually observed on that day. With this sounding as input for the 
model, we expected strong convection to develop. 

Figure 3 shows the characteristics of the modeled storm resulting 
from the simulation at two different simulation times: 15 and 22.5 min. 
The left column contains plots at 15 min, while the right column shows 
the same plots at 22.5 min. The plots from top to bottom represent the 
dynamical and microphysical character of the cloud (see plot caption 
for details), radar reflectivity, vertical velocity, total charge 
density, vertical electric field component, and the horizontal electric 
field component. As is evident from Fig. 3, the model simulation of 
this warm based (~16°C), maritime type storm showed a very rapid 
development with rain formation proceeding via the stochastic coales- 
cence process weli before the ice phase became significant. Note the 
nearly vertical development of the storm due to the lack of wind 
shear. 

The early charge distribution and resulting weak electric fields 
were a result of the limited inductive Interactions between raindrops 
and cloud droplets. Once the Ice phase appeared, the charge separa- 
tion processes became more vigorous and, in the course of a few 
minutes, charge densities reached the order of 10' s of nanocoulombs 
per cubic meter, accompanied by electric field strengths in excess of 
400 kV/m, a value near that necessary to produce lightning breakdown. 
The "classic" thunderstorm dipole structure (positive charge above 
negative) was produced in the simulation with a lower positive charge 
region also In evidence, although there were no observations upon 
which to judge the correctness of the model results. 



$ 


JF 




Fig. 2: Wallops Island, VA (WAL) sounding for 00Z, 13 September 1983. 
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ZiSb. 2 : Characteristics of the Wallops case for two different simulation 
times: left column, 15 min; and right column, 22.5 min. Reading from top 
to bottom, the plots represent: 1) storm dynamical and microphysical 
character including airflow streamlines (dashed lines), cloud boundary 
(solid line), and presence of snow (S), graupel (*), rain (•), and cloud 
ice (-) in amounts greater than a threshold value; 2) radar reflectivity 
in dBz (absolute values systematically too high, but structure is 
representative, contour interval 4 dB at 15 min, 5 dB at 22.5 min; 
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Fly. 3 (continued) : 

3) vertical velocity in m sec -1 (updraft-solid, downdraft-dashed, contour 
Interval 3 m sec -1 ); 4) total charge density in C nr 3 (positive-solid, 
negative-dashed, contour intervals 0.9 pC.m" 3 at 15 min, 1 nC nr 3 at 
22.5 min); 5) vertical electric field in V m" 1 at 15 min, 60 kV m” 1 at 
22.5 min); 6) horizontal electric field (as per the vertical field, 
contour intervals 20 V nr 1 at 15 min; 30 kV m"i at 22.5 min). 
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The model run for this Wallops Island case had to be terminated 
after 22.5 minutes of simulated real time due to the continued buildup 
of the electric fields and the resulting prohibitively small time steps 
necessary to maintain computational stability of the ion transport 
equations. A simulation without electrical processes active was run 
much longer in time. The inclusion of the lightning discharge in the 
model is necessary to overcome this early termination problem. 

Although the electrical simulation had to be terminated 
prematurely, the outputs provided information on the electrical 
environment which might be expected during the early stages of 
thunderstorm development. These data in a form similar to Figs. 3 
and 4, including electric fields, particle densities and charges, and 
ion densities at various altitudes and times within the simulated 
thunderstorm, were sent to EMA for use in conjunction with their 3D 
interaction model (Rudolph et al . , 1984). The plots in Fig. 4 repre- 
sent measurements which an instrumented aircraft would record if it 
made a level pass through the model storm at a given altitude. From 
these plots, EMA scientists were able to extract information con- 
sidered representative of the environment in which the F106 might 
have operated. 

Part of their work involved a determination of the shape factors 
due to the aircraft necessary to interpret the electric field mill data 
which is recorded during each F106 flight. Another aspect of their 
research was to do a parameter study of the aircraft under various 
electrical initial conditions using data from the SEM simulation, where 
appropriate. Their results (Rudolph and Perala, 1985), analyzing the 
data from the electric field mills at the time of a lightning strike to 
the aircraft during flight 83-053, showed that the field was predomi- 
nantly vertical and negative in value. From examination of Fig. 4 at 
22.5 min, the model predicts an environment in the neighborhood of 
28,000 ft which is approaching conditions suitable for lightning 
activity (particle charges in excess of 5 nC/m 3 and electric fields 
approaching 300 kV/m) , indicating that the model is capable of 
generating an environment which is similar to that which was observed. 

The investigation of the low altitude strike problem is aided by 
looking at some of the model results. Figure 5 shows contour plots of 
the vertical electric field component (left panel) with solid contours 
indicating an upwardly directed (positive) field and dotted contours 
representing downwardly directed (negative) fields. The interval of 
each contour is 60 kV/m with the zero line indicated. The right panel 
represents the mixing ratio of graupel in the cloud in grams of graupel 
per kilogram of dry air. The contour interval in this panel is 1 g/kg. 
These plots represent the state of these two model variables for the 
Wallops Island simulation at 22.5 minutes after the initiation of the 
run. 


The three horizontal lines drawn through the contour plots 
represent three possible flight altitudes for the F106. The dashed 
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Fig. 4: Examples of measurements which would be recorded by an 
i nstrumented aircraft flying a horizontal path through the modeled 
storm at 8.4 km AGL and 22.5 min simulation time. Top plot: hail 
(graupel) content in g nr 3 , number of graupel particles nr 1 with 
diameters > 600 gm, and charge on the graupel particles in nC m -3 
(positive charge-solid line, negative charge-dashed line). Bottom 
plot: potential above ground in kV, vertical electric field com- 
ponent (E z ) in dV m" 1 (positive-solid, negative-dashed), and 
horizontal field component (E x ) in kV nr 1 . 
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line at 8.4 km represents the altitude at which the aircraft was 
operating at the time of the lightning strikes. The dash-dot line 
represents a flight level of approximately 10 kft, and the solid line, 
a level of around 19 kft. (The major tick marks on the horizontal and 
vertical axes represent distance in kilometers). These two lower 
altitudes are in the region considered by program personnel to 
constitute the low altitude lightning environment (<20 kft). 

Examination of the figure reveals that, at the altitude at which 
the aircraft was operating on that day, a high electric field 
environment would have been encountered (if the model is adequately 
simulating the cloud which was penetrated). The model results further 
indicate that had the F106 been operating in the vicinity of 10 kft, it 
would have encountered relatively low field strengths and probably not 
have been able to trigger a discharge (the triggering process seems to 
predominate at high altitudes). On the other hand, the flight level 
around 19 kft shows that a very concentrated and intense region of 
electric field in excess of 400 kV/m might have been encountered by the 
F106 had it been flying at that altitude. This result seems to indi- 
cate that lower level, high field regions could exist in the types of 
clouds penetrated during the project; however, examination of the 
right panel in the figure shows that this high field region is coin- 
cident with the maximum hail content of the storm (mixing ratio in 
excess of 14 g/kg) . fn addition, examination of Fig. 3 shows that 
this high field region also corresponds to a region of strong updraft 
and a sharp updraft/downdraft boundary, indicating the possibility of 
strong turbulence. Therefore, while a volume favorable to the ini- 
tiation of lightning by the F106 in the low altitude region may be 
present, the simulation of this one case indicates that the aircraft 
would have to purposely avoid such a region for safety considerations. 
It is also interesting to note that the majority of the low altitude 
strikes that have occurred during the 1985 field season have been in 
the 15-20 kft region. This leads one to speculate that the low alti- 
tude, high field region may be a reasonably frequent feature of the 
storms encountered by the F106 and is not always associated with 
adverse flying conditions. More simulations of the electrical evolu- 
tion of such clouds during their mature and dissipating stages will 
help to clarify this question. 

These results brought to light several questions and limitations 
which resulted in additional research. First of all, while the model 
seemed to be producing realistic simulations of the storm conditions 
within which the F106 was operating and the data provided to EMA was 
useful in their analyses, we had no way of knowing whether or not the 
simulations were actually correct. The only way to verify the ability 
of a model to simulate nature is to have actual observations against 
which to compare the model results. The comparison between the F106 
field mill data and the model output done by EMA scientists was 
encouraging, but was not a sufficient test. Such a test has been done 
on the model for a cold-based continental type storm (Helsdon et al . , 
1984) giving us confidence that the model can correctly simulate cloud 



development and its attendant electrification. However, there are 
significant differences between cold-based, continental clouds and the 
warm based, maritime type clouds found along the east coast. These 
differences are such that we do not feel it is safe to assume that, 
because one type of cloud is well modeled, a different type of cloud 
will be equally well modeled. 

Another, and even more important consideration, involves the 
inability of the model to simulate the electrical development of a 
cloud beyond its early stages. As noted above, the Wallops Island 
simulation had to be terminated after 22.5 minutes simulation time 
because of the constant buildup of the electric field strength. There 
is no mechanism within the model to relax the ever increasing elec- 
trical stresses. Nature takes care of the problem by initiating a 
lightning discharge. In order to be consistent, we realized that it 
was necessary for the model to be able to simulate a lightning dis- 
charge and its effects on the local charges and fields. Only when 
this capability was in place would the model be able to provide simu- 
lations of the mature and dissipating stages of a storm's electrical 
life. Since It is probably rare for the F106 to encounter a storm in 
its youngest stage, the ability to simulate the latter stages becomes 
paramount to further investigations. In order to conduct any further 
investigations into the storm electrical environment, which is a 
desired goal of the project, the model must be capable of simulating 
that environment throughout the life cycle of a storm. 

4.2 Lightning Parameterization 

These considerations led to investigating the question of 
incorporating (parameterizing) the lightning discharge in the SEM. We 
have two types of lightning to deal with; cloud-to-ground and intra- 
cloud. While intracloud lightning is the less well studied, it is 
also the most frequent to occur during a storm. It also seems, 
intuitively, to be the less complicated of the two to incorporate 
in the model and therefore was chosen to be examined first. 

Without going into the arguments pro and con, we chose to adopt 
the philosophy of Kasemir (1960, 1984) whereby the overall Intracloud 
lightning channel is electrically neutral. The charges that exist on 
the channel are the result of ionization along the channel and are not 
"gathered" from the surrounding cloud volume. When one deals with the 
simulation of the lightning discharge, which is a subgrid scale process 
in the SEM, four basic criteria must be established in order to account 
for the process: 1) initiation; 2) direction of propagation; 

3) termination; and 4) charge redistribution. In attempting a first 
order approximation to the lightning discharge, the idea is to keep 
the processes involved as simple as possible while still trying to 
maintain a physically reasonable approach to the problem. Simplicity 
is also in order because the physical mechanisms involved in deter- 
mining these criteria are poorly understood at this point. Thus, we 
have basically employed a single parameter approach to the first three 
processes mentioned above. 



Recently Williams et al t (1985) have made a study of electrical 
discharges in polymethlymethacrylate (PMMA) which had been injected 
with electrons to form various space charge configurations. They found 
adequate scaled correlations to be able to relate the observed behavior 
of the discharges in the PMMA blocks to lightning discharges in thunder- 
clouds. The paper focused on the local electric field strength and 
the space charge distribution as the factors controlling the extent 
and direction of propagation of the discharge channel. While the 
extent and magnitude of the space charge cloud are important in the 
morphology of the laboratory discharges, they identify the local 
electrostatic field as the single most important parameter in 
controlling such discharges. And although there is no direct justi- 
fication for extrapolating their results to what takes place in 
thunderclouds, a scaling approach and the application of their model 
to specific thunderstorm situations argues in favor of accepting a 
possible correlation. As a result of their work and the concept of 
maintaining simplicity in the initial approach to the lightning 
parameterization, it was decided to use the local electric field as 
the parameter controlling the first three criteria in the above list. 

As an initiation criterion, a threshold electric field was chosen 
with a value of 400 kV/m. Williams et al . (1985) quote a value of 
500-1000 kV/m for a breakdown field in thei r Table 1. It has fre- 
quently been mentioned that the in-cloud breakdown field must be in 
the vicinity of 400 kV/m since the highest reliable in-cloud field 
measurements have peaked near that value. We feel that using this 
value in the model is reasonable because the value at a grid point 
represents the average value of the quantity within a grid box 200 m 
on a side (for this model) and this average would most likely consist 
of both higher and lower values. 

The termination criterion was selected to be a critical field 
strength of 150 kV/m based on work done by Griffiths and Phelps (1976) 
involving the propagation of positive streamers in the laboratory. 

While serious questions can be raised about the appropriateness of 
extrapolating such a laboratory value to the atmosphere, we feel that 
it represents a reasonable threshold to use in the first approximation . 
It is true that the stepped leader which precedes a cloud-to-ground 
discharge appears to propagate through regions where the field strength 
does not exceed a few kilovolts per meter and such an extinction 
criterion as proposed above would not seem to be applicable; however, 
we are only dealing with the intracloud discharge at present. 

In keeping with our single parameter formulation for the discharge 
process, we have chosen to use the electric field vector in determining 
the direction of the propagation path. Because the model domain is 
composed of grid points in a rectangular pattern, the lightning path is 
constrained to move either along the side of a grid box or along its 
diagonal, depending upon the direction of the electric field vector at 
that point. An algorithm has been developed which computes the angle 
of the field vector with respect to the vertical direction and 


determines whether the path should be taken along the diagonal or 
the side. Even though the vector may be (and usually is) canted with 
respect to the vertical, the diagonal path is not chosen unless the 
canting exceeds a value of 22.5°. Because of these geometric con- 
straints created by the use of a rectangular grid, the resulting 
lightning channel will probably be artificially relegated to a more 
vertical propagation path than would normally be observed in nature; 
however, for a first approximation, this seems to be a reasonable 
approach which can be improved upon as experience and increased 
knowledge dictate. For example, we recognize that the electric field 
vector at the tip of the propagating leader is important in deter- 
mining the propagation direction; however, this is beyond the scope 
of a first order approximation and is the subject of a current 
research effort. Finally, since the initiation point of the channel 
is the point of highest electric field strength, the channel propa- 
gation is bi-di rectional from that point, moving both parallel and 
anti-parallel to the field vector and terminating when the field 
strength along each segment falls below the cutoff threshold. 

In order to evaluate the effectiveness of these criteria, a test 
was conducted using model output from a simulation of the 19 July 1981 
cloud which was observed in the vicinity of Miles City, MT, and pro- 
duced a single intracloud lightning discharge (inferred from the 
electric field measurements of the NCAR/NOAA sailplane, see Dye et 
al . , 1986). The initiation, propagation, and termination criteria as 
outlined above were applied to the model predicted state of the cloud 
at a time when the electric field was sufficient at some point in the 
domain to exceed the Initiation threshold. From this point, the 
discharge path was computed in two directions following the electric 
field vector. The upward and downward paths continued until their 
termination thresholds were reached, at which point the discharge in 
that direction was stopped. The results of this process are shown 
in Fig. 6 which depicts the model calculated lightning channel 
superimposed on the ambient charge distribution (left panel) and the 
vertical electric field component (right panel). We can see that 
the discharge originates in the region of maximum vertical electric 
field which corresponds to a region of low net charge density. The 
discharge path is basically vertical following the dominant vertical 
electric field component and penetrates both the positive (upper) and 
negative (lower) charge centers. One should note, however, that the 
effect of the horizontal component of the electric field vector (not 
shown) becomes sufficient in the vicinity of the downward propagating 
part of the path so that the channel is deflected to the left as 
termination takes place. The total path length of the discharge 
channel in this instance is 2682 meters. 

The parameterization of the lightning process must be completed by 
specifying the manner in which the charge created along the channel is 
redistributed into the environment. This aspect of the problem is 
being treated using ideas from Kasemir (1984) in which the charge per 
unit length along the channel is proportional to the difference between 
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the potential of the channel (defined as the potential at the point 
where the channel is initiated) and the ambient potential created by 
the cloud charge distribution. Such a scheme guarantees charge 
neutrality of the overall channel (if properly employed) and defines a 
systematic means by which charge can be distributed along the channel 
length. Since the total length of the channel is known, the total 
charge on each segment can be determined. Since intracloud lightning 
may be characterized by filamentary branching into the space sur- 
rounding the main channel, a scheme is also being devised to distri- 
bute the charge created along the channel in the space around the 
channel. This is also required in the model since the numerical 
schemes employed in solving the charge conservation equations are not 
capable of handling the shock of a line discontinuity which would be 
created if the charge were allowed to reside solely along the simulated 
channel. This distribution Is being approximated using a Gaussian 
profile (AexpC-k(x-x 0 ) 2 ]) to distribute the charge horizontally around 
the main channel. Here A represents the magnitude of the line charge 
density along the channel, x 0 is the x coordinate of the channel and k 
is a scaling coefficient chosen to cause the charge value to decrease 
to a specified fraction of its channel value at the edge of the 
distribution region. Since the minimum resolvable scale of a finite 
difference model of this type is twice the grid interval, it is neces- 
sary that a minimum of three grid points be involved in the definition 
of a distribution. In the case of the lightning charge distribution, 
we are starting with four grid points on either side of the channel 
path (1.6 km total width) to establish the distribution region. 
Experience will dictate modifications to this aspect of the problem 
which is the current focus of our research efforts. 



5. CONCLUSIONS 


The results of these investigations have given us some clues to 
the problem of the low altitude strike environment relevant to the F106 
program. In addition, this work has raised several points that have 
lead to additional research thrusts of a pioneering nature. 

With regard to the low level strike environment, the simulation of 
the Wallops Island case has shown that, within the limitations of the 
model physics, a region of high electric field strength does occur in 
the 15-20 kft altitude range. In the early developmental stage of the 
storm, this high field region is associated with the presence of 
graupel , and possibly turbulence, making it a region where flight 
operations would have to be restricted. However, it is entirely 
possible that at a later time in the storm's life cycle, such a high 
field region would continue to exist without being accompanied by 
threatening conditions. In addition, the Wallops Island simulation 
was characterized by little wind shear, leading to a nearly vertical 
development of the modeled storm. In a case where a stronger 
environmental shear were present, the high field region might not be 
concomitant with the presence of graupel or the graupel production 
might not be as heavy as in this case. Such speculation can only be 
resolved by testing the Wallops Island simulation further in time and 
running additional simulations of storms in which environmental shear 
is more significant. 

The ability of the SEM to run a simulation further in time than is 
currently possible depends on the development of the lightning dis- 
charge parameterization. This is the pioneering research referred to 
above. No previous attempts to model the character of the actual 
lightning discharge channel have been made in the context of a large 
scale simulation model of this type. A first order approach to the 
problem has been undertaken in which the ambient electric field is the 
parameter used to determine the initiation, propagation, and termina- 
tion criteria for the simulated discharge. Using these criteria, we 
have been able to create a lightning channel within the context of an 
actual storm simulation. This channel, while not tortuous in nature, 
does exhibit some of the characteristics we would expect to see in an 
intracloud discharge. 

The key to the continued development and use of the SEM as a 
research tool lies in the completion and improvement of the lightning 
parameterization scheme. This work involves the parameterization of 
the charge rearrangement associated with the lightning process, the 
development of a scheme to account for the cloud-to-ground discharge, 
and seeking improvements in the propagation and termination criteria. 
These three aspects of the parameteri zation scheme are currently being 
i nvestigated. 
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